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Transient Convection Due to Imposed Heat Flux: 

Application to Liquid-Acquisition Devices 

Walter M.B. Duval, David J. Chato, and Michael P. Doherty 
National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 

Summary 

A model problem is considered that addresses the effect of heat load from an ambient laboratory 
environment on the temperature rise of liquid nitrogen inside an enclosure. This model has applications to 
liquid acquisition devices inside the cryogenic storage tanks used to transport vapor-free propellant to the 
main engine. We show that heat loads from Q = 0.001 to 10 W, with corresponding Rayleigh numbers 
from Ra = 10 9 to 1 0 13 , yield a range of unsteady convective states and temperature rise in the liquid. The 
results show that Q = 1 to 10 W (Ra = 10 12 to 10 13 ) yield temperature distributions along the enclosure 
height that are similar in trend to experimental measurements. Unsteady convection, which shows self- 
similarity in its planforms, is predicted for the range of heat-load conditions. The onset of convection 
occurs from a free-convection-dominated base flow that becomes unstable against convective instability 
generated at the bottom of the enclosure while the top of the enclosure is convectively stable. A number 
of modes are generated with small-scale thermals at the bottom of the enclosure in which the flow self- 
organizes into two symmetric modes prior to the onset of the propagation of the instability. These 
symmetric vertical modes transition to asymmetric modes that propagate as a traveling-wave -type motion 
of convective modes and are representative of the asymptotic convective state of the flow field. Intense 
vorticity production is created in the core of the flow field due to the fact that there is shear instability 
between the vertical and horizontal modes. For the higher Rayleigh numbers, 10 12 to 10 13 , there is a 
transition from a stationary to a nonstationary response time signal of the flow and temperature fields with 
a mean value that increases with time over various time bands and regions of the enclosure. 

1.0 Introduction 

The investigation of the effect of heat flux applied at the bottom boundary of an enclosure, which 
contains cryogenic liquid nitrogen, to determine its temperature response requires knowledge of the 
dynamical state of the system when there is no applied heat flux; the enclosure is located inside an ambient 
laboratory environment. The no-applied-heat-flux case is necessary since it represents the effect of the heat 
flow from the background laboratory environment on the thermodynamic state of the liquid inside the 
enclosure. As part of a broader program to investigate heat entrapment effects in liquid acquisition devices, 
Bolshinskiy et al. (Refs. 1 and 2) conducted an experiment to measure the effect of background heat flow 
from a laboratory environment on the temperature rise of liquid nitrogen inside an enclosed container 
(Dewar). The results showed that, when no heat flux was applied directly at the bottom boundary, tempera- 
ture rise occurred at various positions along the vertical height of the enclosure. A model problem is 
addressed to obtain insight into the meaning of the temperature measurements (in no-screen, no-heater 
conditions, Ref. 1) for the sole condition of heat flow to the system from the background laboratory 
environment. 

To determine a suitable range of parameters, we estimated a range of typical heat loads caused by the 
thermal background radiation, which yielded a conservative heat load of Q = 0.001 W to a typical maxi- 
mum value of Q = 10 W. We approximated the problem as a rectangular cross section of a cylinder; this 
was convenient since, at the lowest heat load (Q = 0.001 W), the Rayleigh number Ra is on the order of 
10 9 and flow asymmetry occurs in this system for Ra on the order of 1 0 7 , as shown by Duval, Chato, and 
Doherty (Ref. 3). This implies a full transient calculation for a three-dimensional cylinder. However, a 
rectangular model provides an approximation to the system that allows the transient dynamics to be 
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investigated beyond asymmetry. The dynamics of the system for Ra = 1 to 1 0 8 has been considered for an 
aspect ratio of Ar= 1 (Ref. 3), and a short-time solution for Ra = 10 9 is presented for Ar = 3.24 (Ref. 3). 
Herein, we consider the asymptotic dynamics for Ra = 10 9 and Ar = 3.24, and we extend the calculations 
to Ra = 10 13 for the highest heat-load condition. We show that at the highest heat load, the model predicts 
temperature magnitude comparable to that in the experiments (Refs. 1 and 2). In addition, the model 
shows the dynamical state of the accompanying flow field, which was not accessible during the 
experiments since only temperatures inside the system were measured (Refs. 1 and 2). 

In analogy to our model problem of a cryogenic liquid inside an enclosure that has a large height and 
that absorbs heat flux predominantly from its bottom and top boundaries, there have been a number of 
investigations on the convective instability inside enclosures heated from below with prescribed temperature 
boundary conditions on the bottom and top boundaries, which lend insight into our model problem. 
Townsend (Ref. 4) showed that convection near a heated horizontal surface depends on the layer near the 
heated surface (independent of the geometrical height of the enclosure) and is determined by the heat flux, 
viscosity, and conductivity of the fluid. Townsend described the formation of thermals on the heated surface 
as rising columns of hot air from the edge of the conduction layer. A phenomenological theory was put forth 
by Howard (Ref. 5) who postulated that convection consists of a conductive phase, during which the 
conduction boundary layer is formed, followed by a short time interval, during which the thermals break off. 
Convective instability, or the formation of thermals, occurs after a time f that is shorter than the thermal 
diffusion time H 2 / a, where H is the enclosure height and a is the thermal diffusivity. However, the thermals 
break away at a time that is short in comparison to the conductive period f. In addition, Howard (Ref. 5) 
showed that the mean temperature profile, based on the idea of marginal stability of the mean flow, depends 
on the thermal boundary layer thickness 8 when it becomes unstable (in which the local Rayleigh number is 
based on 8) and shows reasonable agreement with experimental data. 

A number of two-dimensional numerical and experimental investigations of convective instability of 
flows inside enclosures that are heated from below consider the short-time evolution of the instability (Refs. 6 
to 9) instead of the statistically steady regime (Ref. 1 0) for long times. Elder (Ref. 6) showed that a Hele- 
Shaw cell can be used to simulate, experimentally and computationally, two-dimensional flow in a porous 
medium within the context of flow in a homogeneous horizontal slab that is heated from below. When the 
approximation of a Hele-Shaw cell was used, numerical prediction of unsteady convection (Ref. 7) inside an 
enclosure in which the width of the heating zone at the base can be adjusted showed that the first motion is a 
set of eddies growing on opposite ends of the local heating zones. As the instability evolves, additional eddies 
are formed in the central region of the heating zone. For a long heating zone that is approximately the hori- 
zontal length of the slab, the initial motion of a pair of end cells followed by successive growth of a string of 
cells above the heater is similar to the case of a short heating zone. These events are phenomenologically 
similar to the short-time prediction of our model for the initiation of convective instability. 

For a gravitationally unstable layer of viscous fluid initially hot below and cold above a horizontal 
plane, Elder (Ref. 8) showed that the scale of the motion is set by the thickness of the interface 8 and is 
independent of the vertical distance H of the horizontal plane in agreement with Townsend’s (Ref. 4) 
findings. The amplification process of the disturbance in the interface is shown to be similar for the flow 
through a porous medium and for a viscous fluid inside a Hele-Shaw cell with impermeable boundaries. 
Experimental observation (Ref. 8) confirming numerical prediction shows that the instability develops as an 
array of blobs above a lower surface in a proto-sublayer. The disturbance remains embedded in the proto- 
sublayer until it reaches a finite amplitude. Similar events for short-time dynamics are shown to occur in our 
model for Rayleigh numbers on the order of 10 12 . The effect of various types of disturbance heat sources at 
the bottom boundary — namely, random, spatially varying, and constant — show similar developments of the 
instability (Ref. 8). This implies that the instability can be induced with a constant heat source as shown in 
our computational model. The development of the proto-sublayer is shown to be similar to the phenomeno- 
logical description of the sublayer proposed by Howard (Ref. 5): that the mean temperature profile depends 
on the boundary layer thickness of the sublayer in which conduction effects dominate. 

Elder (Ref. 9) shows that using one-dimensional mean field equations to model the flow development of 
the proto-sublayer can simulate two-dimensional numerical models. The spatial means are taken over a 
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horizontal plane in the mean field equations, and the simulation is shown for Rayleigh numbers on the order 
of 10 5 . Even though Elder’s numerical simulation (Refs. 7 and 8) captures the growth of cells in the flow 
field as the convective instability evolves, the temperature field shows that thermals are roughly represented 
by blobs; this is also shown experimentally for two-dimensional flows simulated with a Hele-Shaw cell 
configuration. The resolution of mushroomlike thermals exhibiting a blunted nearly hemispherical cap was 
captured in the experiments by Sparrow, Husar, and Goldstein (Ref. 1 1), which used water inside a three- 
dimensional enclosure heated along a horizontal surface. They showed that thermals are generated 
periodically in time at fixed sites and that the breakup Rayleigh number based on the boundary layer 
thickness of the conduction layer is a constant, which validates a prediction by Howard (Ref. 5). Even 
though, for the short-time development of convective instability (in particular the proto-sublayer (Ref. 8)), 
the preceding works show trends that are similar to those of our model, for the long-time-scale asymptotic 
dynamics, we show herein that the dynamical state of the system is a traveling-wave-type motion of 
convective modes, which damps the evolution of thermals near the bottom boundary. 

In the following sections, the formulation of the basic physics of the problem and discussion of the 
numerical solution are presented. Then, the global dynamics of the flow field — showing the growth of 
thermals, the transition to asymmetry, and the traveling-wave -type motion of the convective modes — are 
presented. Next, the flow field is probed at various fixed points to show time histories of the flow and 
temperature fields. Finally, we identify the range of heat load from the ambient laboratory environment, 
showing predicted temperature magnitudes that are comparable to those of the experiments. 

2.0 Formulation 

The model, which addresses the effect of the external heat load Q from an ambient laboratory 
environment at temperature T x = 293 K on the temperature rise of saturated liquid nitrogen, initially at 
T a = 77 K, inside an enclosure is shown in Figure 1. The heat load transmits constant heat flux q" on 
the boundaries of the enclosure with equal heat flux absoiption on the bottom and top ( q " = q"), whereas 
the sidewall heat flux q" is assumed to be small, q" « q h ", since insulation was used at the side 
boundaries in these experiments (Refs. 1 and 2) to investigate heat entrapment effects inside liquid 
acquisition devices. The figure shows the locations of points inside the enclosure ( x p , y p ) that allow 


Q qf Q 



Figure 1 . — Physical description of cryogenic fluid (ambient temperature, Ta, 
and density, p^) inside an enclosure absorbing surface heat flux, at the 
bottom, sides, and top ( qt, ", q s ", and q t "), from ambient environment heat 
flow, Q. Location of points in the flow field are denoted as ( x p , y p ). 
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probing of the time history of the flow and temperature fields similar to the experimental measurement 
points. The maximum heat load absorbed from the ambient owing to background thermal radiation can be 
estimated to be on the order of 10 W, whereas a minimum conservative value is assumed to be 0.001 W. 
The absorption of an infinitesimal heat load Q through the boundaries of the large enclosure with height 
H= 61.9 cm, approximately 3 times its horizontal length, L= 19.1 cm, drives the system off equilibrium. 
A deviation from equilibrium conditions will cause an increase in temperature given by 

DT 

= aV 2 T (1) 

Dt 

which can give rise to fluid motion owing to its coupling with the body force (g v = ng 0 ) through the 
density field. The body force is directed downward vertically, in which n represents the various 
gravitational levels of the Earth (n =1), Moon (n = 1/6), or Mars (n = 1/3), and g G is the gravitational 
acceleration on Earth. For an incompressible Boussinesq fluid (where V is the velocity vector field), 

V • V = 0 (2) 

When the density p is independent of pressure but is a function of temperature, the density can be given 
as an equation of state: 


p = p(T 4 )[l-pA7] (3) 

where P = -(l/p)(op/f 7) /; is the coefficient of thermal expansivity, the subscript p indicates pressure, and 
A T= T - T a is the temperature difference. The motion of the fluid satisfies 

DV Vo 

~~r— = — + vV 2 V - pAT «g 0 (4) 

Dt p 


where t is time and p is pressure. 

The dynamical equations of motion (Eqs. (1), (2), and (4)) satisfy the initial condition that liquid 
nitrogen is at saturation: 


t = 0, T(pc, y, 0) = T a (5) 

as well as the no-slip boundary condition along the boundary T of the enclosure, 

t> 0, V = 0 on T (6) 


Prescribed heat flux conditions are imposed along the bottom and top boundaries T of the enclosure: 


y 


= 0, Vx, q" = ql, dT / 


y - H,\/x, q" = q;, 
and along the left and right walls, 


x - 0, Vy, q’ = q', 


/ dy 


- _qb / 


y = o 


K 


dT, 


/ dy 




y=H 


K 


dT, 


/ dx 


= - V 

x=() /K 


( 7 ) 

( 8 ) 

( 9 ) 
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, = l,vv, 


( 10 ) 


where k is the thermal conductivity. 

Since the top heat flux q" has been shown to stabilize the flow (Ref. 3), fluid motion is generated by 
either natural convection due to the wall heat flux q" or by convective instability due to the bottom heat 
flux ql . Since there is a threshold value on convective instability determined by q” b prior to its onset, a 
natural convection base flow is always generated owing to the finite heat absorption q" at the wall. Thi s 
base flow is independent of the heat flux absorbed at the top and bottom boundaries. 


3.0 Scaling Analysis and Computational Method 

For a computational solution of the dynamical equations of motion (Eqs. (1), (2), and (4)), the vector 
momentum equation (Eq. (4)) can be transformed to a scalar equation for the vorticity field for two- 
dimensional flow. A single vorticity component in the z direction is defined as 


8v du 
dx dy 


(11) 


where v and u are the vertical and horizontal components of velocity, respectively. The resulting vorticity 
equation — when scaled using the characteristic length, time, and velocity (L, L 2 /a, and a/L ) — becomes 


d^_ 

dt* 


+ 


u — + v — 
dx dy' 


—Pr 


dX d 2 d, 




dx* 2 dy 


*2 


+ Ra Pr 


J 


dT* 

dx 


( 12 ) 


which implies that vorticity is generated by the horizontal temperature gradient that causes its advection 
and diffusion. Elere, Pr is the Prandtl number and Ra is the Rayleigh number; a superscript asterisk (*) is 
used to denote a dimensionless quantity. The continuity equation (Eq. (2)) is satisfied identically when the 
velocity components are expressed in terms of the gradient of the scalar stream function \\i in which 
u = d\\i/dy and v = - d\\i/dx; this transformation, when applied to the vorticity component equation 
(Eq. (11)), yields Poisson’s equation: 


d 2 \\i* d 2 \\i* 

dx* 2 dy* 2 


(13) 


which is used to solve for the flow field (u, v). The flow field is generated by the buoyancy source term in 
the vorticity field equation from the diffusion of heat along the boundary of the enclosure communicated 
through the density field. The buoyancy source term is coupled to the temperature field given by its 
dimensionless form: 


dT* f * dT* t dT*) 

11 TT^ + V 

dt ^ dx dy J 


' d 2 T* 


d 2 T* s 
~dy^^ 


obtained using the following temperature scale: 


t-t a 

q" L/k 


(14) 


(15) 
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In addition to Ra and Pr in the vorticity equation (Eq. (12)), which are defined as 


Ra= MmoV 

Kva 


and 


Pr = — 


a 


(16a,b) 


(where v is the kinematic viscosity), the remaining parameters of the problem are obtained by nondimen- 
sionalizing the boundary conditions. This produces the heat flux ratios for the bottom and top boundaries, 
fb and /: 


y 


= o, dr 


r * 
oy 


l/=o 


= ~f b 


and 


y* = Ar, dT 


<-\ * 
oy 


=-f, 


I y*=Ar 


and the heat flux ratios of the sidewall boundaries,/,: 


= o, sr 


* 

OX 


lx*=0 


=-/; 


and 


= 1 8T* 


^ * 
ox 


U =1 


=-/; 


(17) 


(18) 


(19) 


( 20 ) 


An additional four parameters are obtained. These are the aspect ratio Ar and dimensionless heat fluxes f b , 
f, and / defined as 


, H 
Ar — — 
L 

(21a) 

rr 

f Qb 
Jb ~ — 

(21b) 

< l 


f q” t 

(21c) 

q 


Js „ 

(2 Id) 


The heat flux scale /'used in Equations (21b), (21c), and (2 1 d) depends on the dominant prescribed 
heat flux at the enclosure boundary that causes convection. For this problem, it is the bottom heat flux 
that drives convective instability since q" « q". In general, the parametric set A consists of six 
parameters: 


A = A (Ra,Ar,Pr,f,f b ,f) 


( 22 ) 


To simplify the parametric set, we consider a specific cryogenic fluid with fixed geometry such that 
Pr and Ar are constant. In addition, since we assume very low heat leakage at the sidewalls,/ is also 
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constant. The heat load from the ambient environment affects Ra,f b , and /. Since f b =/ and the use of the 
heat flux scale q" implies that /, = 1 and f = 1 , variation of the heat flux appears in/,. Thus, the reduced 
parametric set A' becomes 


A' = A \Ra,f) (23) 

This implies that the effect of heat load from the ambient laboratory environment (ng 0 = lg 0 ) on 
the temperature rise of liquid nitrogen inside an enclosure can be characterized by varying only two 
parameters: the Rayleigh number and the sidewall heat flux ratio. This simplification yields a reasonable 
set of parameters for computational parametric study. 

To obtain a computational solution of the dynamical equations of motion for the reduced parametric 
set, we employ finite difference methods. First, we use a direct solver approach to solve Poisson’s 
equation (Eq. (13)) using a matrix inversion that yields the flow field. The input of the flow field is used 
to solve for the temperature field in Equation (14) via a flux-corrected transport method used to resolve 
small-scale thermals. Finally, since the horizontal temperature gradient drives vorticity production 
through buoyancy, we solve the vorticity equation (Eq. (12)). 

This direct solver approach requires the solution of Poisson’s equation for the first time step only. 
Subsequent time steps use the vorticity field (Eq. (12)) as input to the temperature field to obtain a time- 
marching solution of the flow and temperature fields. Since we are interested in time-accurate solutions, 
we use a third-order Adams-Bashforth time-discretization scheme with a typical time step of approxi- 
mately one -thousandth of a second: At = 0.00125 s. Such a small time step allows an asymptotic solution 
of the flow field to be computed. Since we are considering relatively high Rayleigh numbers, 10 9 to 10 13 , 
the smallest grid size employed is 1000 by 1000. Owing to the large vertical height of the enclosure, 
asymptotic dynamics of the flow field require long-time computations, which can take on the order of 500 
to 1000 hr on a Silicon Graphics International Corporation (SGI) Altix supercomputer system. Memory 
requirements for each case range from a low of 12 CPU to a midrange of 40 CPU and a high of 100 CPU 
on supercomputers with capacities of 508, 1024, and 2048 CPUs (1 CPU is approximately 2 GB). Thus, 
substantial computational resources are required to obtain solutions for the reduced parametric set. 

4.0 Numerical Results 

4.1 Global Dynamics of the Flow Field for High Rayleigh Numbers 

Owing to the large geometric length scales of the enclosure (L = 19.1 cm and H= 61.9 cm), com- 
putation at the lowest Rayleigh number, for the most conservative estimate of the heat load from the 
background laboratory environment, is on the order of 10 9 . For the range of Rayleigh numbers (10 9 to 
10 13 ) considered to investigate the effect of background heat load Q from 0.001 to 10 W, we show that 
departure from thermodynamic equilibrium ensues from convective instability. The instability propagates 
from the bottom to the top of the enclosure and asymptotes to a traveling-wave-type motion of convective 
modes toward dynamical equilibrium. The global dynamics of the flow field indicate self-similarity of the 
flow as an asymptotic state is approached. Elowever, the approach to the asymptotic state depends on the 
magnitude of the Rayleigh number. 

For the most conservative estimate of background heat load ( Q = 0.001 W, which corresponds to 
Ra = 10 9 ), Figure 2 shows a typical base flow from which convective instability ensues. In response to heat 
load from the ambient, a heat flux of q", - q" = 3.49 x 10 1 erg/cm 2 -s is transmitted through the top and bot- 
tom surfaces, and careful insulation of the sidewalls allows for a minimal heat flux of q" = 0.001 erg/cm 2 -s. 
Throughout the parametric variation (the increase in background heat load), the absoiption of heat flux 
through the sidewall is maintained constant at q" = 0.001 erg/cm 2 -s. Figure 2 shows that prior to convective 
instability, the sidewall heat flux/ causes natural-convection-dominated base flow with two large counter- 
rotating cells. This typical base flow exists for all Rayleigh numbers. The temperature field shows the 
conduction boundary layer 8 built up at the top and bottom of the enclosure owing to heat flux / and /,. This 
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is similar to the conduction sublayer in the phenomenological theory of Howard (Ref. 5) and the 
proto-sublayer of Elder (Ref. 8). 

Figure 3 shows the transition by which the base flow field becomes unstable against convective 
instability. The instability of the thermal boundary layer thickness 8 is evident at t = 300 s and shows that 
the disturbance remains embedded in the proto-sublayer as pointed out by Elder (Ref. 8). Due to the fact that 
heat flux is absorbed at the top and bottom of the enclosure, convective instability unfolds at the bottom of 
the enclosure owing to an adverse temperature gradient stemming from the unstable density gradient driven 
by buoyancy, while the top of the enclosure is convectively stable since the density gradient is stably 
stratified. These results are in agreement with the findings of Townsend (Ref. 4) and Elder (Ref. 8) that 
convective motion at early times is dependent on the thermal boundary layer thickness and is independent of 
the height of the enclosure. As shown in Figure 3, the six cells that are generated initially are segregated (at t 
= 300 s) toward the bottom wall comer, similar to the findings of Elder (Ref. 7). The growth of these cells 
yields small-scale thermals (at t = 432 s) (as shown by the temperature field) while the six cells self-organize 
into two larger cells with small secondary cells (as shown by the flow field). The amalgamation of the cells 
indicates a local increase in the wavelengths of the cells as also observed by Elder (Ref. 7). According to 
Elder (Ref. 8), the amalgamation of cells was also observed in laboratory experiments and represents a 
broadening of the spatial spectrum of the velocity field to lower wavenumbers or higher wavelengths. 

The transition from short-time -scale events, in which thermals break off, to long-time -scale events, 
which mark the propagation of the instability, occurs for t > 432 s. The shear instability of the top 
convective cell, with the quiescent fluid on top, fomis a local source that excites the growth of secondary 
vortices ( t = 432 s). Self-repetition of the internal shear instability mechanism leads to the propagation of the 
local convective instability front toward the top of the enclosure ( t = 1728 s). As the top of the enclosure is 
approached, this system of four symmetric counterrotating vortices becomes asymmetric, and the break- 
down of symmetry in the system is shown for t = 1872 s. Regions of shear instability between the four cells 
bifurcate to a series of secondary cells, which leads to self-organization of the flow field into a set of four 
asymmetric cells. These cells are precursors to a traveling-wave -type motion of convective modes. The 
traveling-wave-type mode occurs after the computation time of 1872 s. However, as will be shown, 
increasing the background heat load Q to 0.01 W or Ra to 10 10 increases the speed at which the instability 
propagates toward the top of the enclosure. Note that in Figure 3, and in the following figures, the local 
coordinates, x p and y p , are used for graphics which are normalized, for convenience, with respect to L and H, 
respectively. However, the normal dimensionless coordinates, x and y , which bound the internal region of 
the enclosure, are normalized with respect to the horizontal length scale L. 
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Figure 2. — Base flow field, due to sidewall heat flux, which becomes unstable owing to local 
convective instability at the bottom boundary, where vy is the stream function. Heat flux ratios at the 
bottom, top, and sidewall, f b - 1 , f t = 1 , and f s = 2.87x1 0 -5 , respectively; aspect ratio, Ar = 3.24; 
Prandtl number, Pr= 2.27; Raleigh number, Ra = 10 9 . 
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Figure 3. — Global dynamics of the flow field, showing the growth of thermals over time, f, and the 
propagation of the flow field leading to incipient asymmetry; heat fluxes at the sidewall, bottom, 
and top, g" = 0.001 erg/cm 2 -s and q" b =q”= 3.49x1 0 1 erg/cm 2 -s, respectively; heat flux ratios at 

the bottom, top, and sidewalls, fb = 1 , ft - 1 , and f s = 2.87x1 0 -5 , respectively; heat load, Q = 0.001 
W; aspect ratio, Ar= 3.24 (height and width of enclosure, H = 61.9 cm and L = 19.1 cm); Prandtl 
number, Pr= 2.27; and Raleigh number, Ra - 10 9 . Results are shown for increasing time. 
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4.2 Traveling- Wave-Type Mode 

The effect of increasing the absorption of heat on the boundary of the enclosure from the ambient 
environment is to increase the speed of transition to a traveling-wave-type motion of convective modes. 
For the same time scale as in Figure 3 (t= 1872 s), Figure 4 shows for Ra = 10 10 the propagation of the 
basic four symmetric modes ( t = 432 s) and their breakdown to asymmetric modes (t = 576 s), which 
leads to the traveling-wave-type motion of convective modes ( t = 864 to 1872 s). Unlike for Ra = 10 9 , the 
transition to asymmetric modes (t = 576 s) occurs before the instability front reaches the top of the 
enclosure. The asymmetric front t = 576 s (t< 792 s) propagates toward the top of the enclosure through 
shear instability, which generates secondary modes until the top of the enclosure is reached (t = 792 s) 
and a fourth mode emerges. The fourth mode is a precursor to a traveling-wave -type mode ( t = 864 s) that 
propagates between the top and bottom of the enclosure as time increases. As shown by the vorticity field 
in Figure 4, the dynamics of this traveling-wave-type mode generates internal vorticity that mixes the 
liquid. Consequently, the temperature is uniform along the height of the enclosure even though, as will be 
shown, the flow and temperature fields are oscillating in time locally ( t = 1872 s). 

The characteristic observation time for the laboratory experiments for the case of no heater and no 
screen (Ref. 1) is t = 600 s. Figure 5 shows the characteristics of the flow and vorticity fields for a time 
scale of 720 s in order to gauge typical convective events. Dynamical events in Figure 5 are shown as a 
function of increasing heat absoiption Q from 0.001 to 1.0 W, corresponding to increasing Ra from 10 9 to 
10 12 . The effect of increasing the background heat load Q is to increase the speed at which the instability 
propagates toward the top of the enclosure and its transition to a traveling-wave-type mode. Figure 5 
shows that the shear instability between the asymmetric modes generates the production of internal 
vorticity in which the intensity increases with increasing Rayleigh number. As Ra increases from 10 10 to 
1 0 12 , the slight decrease in the maximum temperature near the top of the enclosure is evidence of the 
effectiveness of mixing due to the intense production of vorticity. On the other hand, there is a fast 
increase in the maximum velocity — from 0.83 to 7.48 cm/s. Even though the background heat load from 
the experiment (Ref. 1) is unknown, Figure 5 gives an indication of the effect of a range of heat-load 
conditions on convective instability. For Ra = 10 9 and Q = 0.001 W, the instability penetrates only about 
20 percent of the height of the enclosure over the time scale of the experiment. Thus, it is unlikely that a 
heat load of Q = 0.001 W would be absorbed. The temperature measurements from the experiments 
indicate that convective instability has penetrated the entire height of the enclosure, implying that 
Q> 0.1 W. Even though the depth of penetration of the instability helps to establish a lower bound on the 
heat absoiption, it is the magnitude of the increase in temperature with time that dictates the neighbor- 
hood of the heat-load absoiption. This will be shown by probing the flow field at various fixed points 

( x pi yp)- 


4.3 Small-Scale Flows 

As the background heat load increases beyond Q = 0.01 W (Ra = 1 0 10 ), it becomes challenging to 
obtain computational solutions owing to the increase in mode numbers, which characterizes the onset of 
convective instability. Figure 6 shows the bifurcation to higher mode numbers for increasing heat load or 
Rayleigh number. The first bifurcation — from six to eight modes — occurs for Ra = 10 11 , which is charac- 
terized by intensive vorticity production near the bottom boundary. The second bifurcation — from 8 to 1 0 
modes — occurs for a decade increase in Rayleigh number, Ra = 10 12 , or heat load, Q = 1 W. After another 
decade increase in Rayleigh number, Ra = 10 13 , the mode numbers double to approximately 20 modes. 
This explosive increase in mode numbers required an increase in grid size to 1500 by 1500 for small-scale 
resolution. The explosive increase in mode numbers caused a considerable challenge in obtaining compu- 
tational solutions of time -asymptotic dynamical states for Ra = 10 13 . Small-scale thermals erupt from the 
bottom wall on the same time scale as the growth of the instability; this is shown for Ra = 10 13 and 10 11 in 
Figure 6. 
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Figure 4. — Transition to asymmetry leading to the traveling-wave-type motion of convective modes; 
heat fluxes at the sidewall, bottom, and top, g" = 0.001 erg/cm 2 -s and ql = q " = 3.49x1 0 2 erg/cm 2 -s, 

respectively; heat flux ratios at the bottom, top, and sidewalls, fb = 1 , ft - 1 , f s - 2.87x1 0 -6 , 
respectively; heat load, Q = 0.01 W; aspect ratio, Ar= 3.24 (height and width of enclosure, 

H = 61.9 cm and L = 19.1 cm); Prandtl number, Pr= 2.27; and Rayleigh number, Ra = 10 10 . 

Results are shown for increasing times, t. 
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Figure 5. — Effect of heat flux absorption at the sidewall, bottom, and top ( g" = 0.001 erg/cm 2 -s and 
q" = q” = q". respectively) on vorticity production over a time scale of t = 720 s. Heat flux ratios 

fb= 1, ft - 1, f s = 2.87x1 O' 5 to 2.87x10^ (corresponding to heat loads of Q = 0.001 to 1.0 W); 
aspect ratio, Ar = 3.24 (height and width of enclosure, H = 61.9 cm and L = 19.1 cm); and Prandtl 
number, Pr= 2.27. Results are shown for increasing heat flux absorptions, g"; heat loads, Q; and 
Rayleigh numbers, Ra. 
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Figure 6. — Bifurcation to higher mode numbers for increasing Rayleigh numbers (grid size increases from 
1000 by 1000 to 1500 by 1500 for Rayleigh number, Ra = 10 13 ) and small-scale thermals. Aspect ratio, 
Ar = 3.24; and Prandtl number, Pr = 2.27. Results are shown for different times, f; heat fluxes, g"; heat 
loads, Q; heat flux ratios at sidewalls, f s ; and Rayleigh numbers, Ra. 


Since the resolution of small scales becomes important for the highest Rayleigh number, Ra = 1 0 13 , or 
heat load, Q= 10 W, considered, we show in Figure 7 the effect of grid size, which illustrates sensitivity to 
initial conditions. At Ra = 10 13 , the system is highly nonlinear and, as will be shown, becomes chaotic for 
Ra < 10 13 . A chaotic dynamical system implies sensitivity to initial conditions. For a time scale of t = 9 s, 
the increase in grid numbers from 1000 by 1000 to 2500 by 2500 serves as a test for sensitivity to initial 
conditions. The increase in grid numbers represents a perturbation to the base flow solution and shows slight 
nuances in the dynamics of the modes; however, a grid size of 2000 by 2000 is adequate to resolve the 
short-time-scale dynamics of the growth of modes. The difference between 2000 by 2000 and 2500 by 2500 
is due to the fact that the dynamics of mode interactions is also sensitive to the grid number density. In 
addition, self-organization of the modes occurs rapidly, which can result in slight nuances in mode numbers. 
This implies that a pattern similar to that for the 2000 by 2000 grid occurs for a 2500 by 2500 grid, but at 
t < 9 s. Because of this sensitivity to the initial condition, the exact patterns are not be repeated at fixed time 
t = 9 s as the grid size increases. Flowever, certain basic patterns remain the same: that is, the explosion to 
higher mode numbers and the increase in vorticity production near the bottom wall. 

Figure 8 shows the self-similarity of the propagation of the convective instability for a 2000 by 2000 
grid at Ra = 10 13 , which corresponds to the upper bound in heat load. This figure illustrates that, beyond 
the onset of convective instability at t = 8 s, in which there is an explosion in mode numbers, the growth 
of the instability, 9.6 s <t< 19.2 s, occurs through the self-organization of modes, which results in the 
basic four symmetric counterrotating vortices that also occur for the lowest Rayleigh number (Ra = 1 0 9 ). 
Beyond t > 19.2 s, this convective state follows trends similar to those shown for the lower Rayleigh 
numbers: that is, transition to asymmetry leading to a traveling-wave-type mode. As the flow propagates 
upwards, there is an increase in the small-scale internal production of vorticity in the bulk flow field. 
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Figure 7. — Illustration of sensitivity to initial conditions showing the effect of grid size on the resolution of 
small scales at time, t = 9.0 s, for Rayleigh number, Ra = 10 13 ; heat flux, g" = 3.49x10 s erg/cm 2 -s; heat 
load, Q = 10.0 W; heat flux ratios at the bottom, top, and sidewall, fb - 1, ft = 1, f s = 2.87x1 0~ 9 , 
respectively; aspect ratio, Ar= 3.24; and Prandtl number, Pr= 2.27. 


Stream function, Vorticity, 




Figure 8. — Evolution of instability at small scales, showing the self-organization of modes and vorticity 
production as the instability propagates from time, f, of 8 s to 1 9.2 s. Heat flux, q" = 3.49x1 0® erg/cm 2 -s; 
heat load, Q = 10.0 W; aspect ratio, Ar= 3.24; Prandtl number, Pr= 2.27; Rayleigh number, Ra = 10 13 ; 
and grid size, 2000 by 2000. 
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4.4 Scalar Measures 


For the time scale of the experiments (Ref. 1), Figures 9(a) and (b) show the characteristic maximum 
velocity and temperature due to heat absoiption from the ambient background environment. There is a shaip 
increase in the maximum velocity as the heat absorption or the Rayleigh number increases (Fig. 9(a)). The 
fast increase in velocity is the cause for the decrease in time scale for the propagation of the convective 
instability along the height of the enclosure. An increase in the heat absoiption from the environment causes 
a rise in the maximum temperature at the top boundary (Fig. 9(b)). For a short time (t = 360 s), the maxi- 
mum temperature remains nearly constant, and for a long time (t = 720 s), the maximum temperature 
decreases as the Rayleigh number increases due to the fact that intensive convective motion mixes the fluid 
in the top boundary layer near the wall. 



Rayleigh number, Ra 

Figure 9. — Maximum velocity (l/ ma x) and temperature (F ma x) as a 
function of Rayleigh number, Ra. Aspect ratio, Ar- 3.24; Prandtl 
number, Pr= 2.27; heat flux ratios at the bottom and top, f b = 1 and 
ft - 1 , respectively, (a) Maximum velocity in flow field as heat 
absorption from the ambient environment increases for a fixed time, 
(b) Effect of heat absorption from the ambient environment on 
maximum temperature inside the enclosure for a fixed time. 
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Nondimensionalized temperature, T 

Figure 10. — Temperature profile along the vertical axis atx p = 2.87 cm and 
fixed time t = 648 s for increasing Rayleigh numbers at ambient conditions 
for two heat loads: heat flux ratios at the sidewall, f s = 2.87x1 CT 5 (heat flux, 
q" = 3.48X10 1 erg/cm 2 -s; heat load, Q = 0.001 W) and f s = 2.87x10" 8 
(g" = 3.49x1 0 4 erg/cm 2 -s, Q = 1.0 W). Aspect ratio, Ar= 3.24; Prandtl 
number, Pr= 2.27; and heat flux ratio at the bottom and top, fb = 1 and 
ft =1, respectively. 

The temperature profile along the vertical axis is contrasted for low Q = 0.001 W (Ra = 1 0 9 ) to mod- 
erate Q = 1 W (Ra = 10 12 ) heat absoiption in Figure 10. Over the time scale of the experiment (Ref. 1) 

(t = 600 s), the propagation of convective instability penetrates about 20 percent of the height of the enclo- 
sure for Q = 0.001 W, whereas the instability propagates the entire height of the enclosure for Q = 1 W. For 
Ra = 10, the temperature remains nearly constant throughout the bulk of the liquid, with a large tempera- 
ture gradient in the stably stratified layer near the top wall and a well mixed layer near the bottom wall. 


4.5 Local Dynamics at Fixed Points 


The response of the temperature and flow fields at a fixed point (x* = 0. 1 5, y* p = 0. 1 ) near the bottom 

of the enclosure for long (Ra = 10 9 to 10 12 ) and short (Ra = 10 13 ) time bands for various heat loads from the 
ambient laboratory environment is shown in Figure 11. These predictions indicate the response time for the 
onset of convective instability as the Rayleigh number or heat absoiption increases. The response time of 
the instability at 1 0 percent of the enclosure height decreases from t = 400 s for Ra = 1 0 9 to t = 15 s for 
Ra = 10 13 as more heat is absorbed from the environment. Both the temperature and velocity fields oscillate 
aperiodically from Ra = 10 9 to 10 13 , this is due in part to the large dimensions of the enclosure, Ar = 3.24 
(H = 61.9 cm and L = 19.1 cm). The velocity and temperature fluctuations increase as the heat load or 
Rayleigh number increase. The time series for the velocity field component is nearly stationary in the mean, 
varies about a fixed level for Ra = 1 0 9 , and varies slightly about the mean for a segment of the time series 
for Ra = 10 10 and 10 11 . However, for Ra = 10 12 to 10 , the velocity field is nonstationary in the mean, does 
not vary about a mean fixed level, and has an upward trend. This implies that the mean varies with time for 
high Rayleigh numbers. An estimation of an average value from such a time series should take into account 
the variation of the mean with time. On the other hand, the temperature field is stationary for Ra = 10 9 to 
10 10 , and a nonstationary time series occurs for Ra = 10 11 and 1 0 1 2 for t > 200 s, which started out as 
stationary for the short time band t < 200 s. Since we show only a short time band of the time series for 
Ra = 10, stationary behavior appears; however, a nonstationary trend would ensue for a time band on the 
order of 800 s, similar to Ra = 10 . These trends indicate that the averaging of temperatures for various 
points on a plane should take into account the nonstationary behavior of the time series so that the 
temperature histories can be correlated properly with the convective states of the flow field. 
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Figure 1 1 . — Time history of nondimensional temperature, T , and velocity fields (horizontal component of 
velocity, u) near the bottom of the enclosure (x* =0.15, y* =0.1) , showing an increase in intensity of the 
fluctuations as the Rayleigh number, Ra (heat absorption) increases for aspect ratio, Ar= 3.24 (height and 
length of enclosure, H = 69.1 cm and L = 19.1 cm); Prandtl number, Pr= 2.27; and heat flux ratios at 
bottom and top of enclosure, ft, - 1 and f t = 1 , respectively. Results are shown for different heat flux ratios 
at the side wall, f s , heat fluxes, g", heat loads, Q, and Rayleigh numbers, Ra. 
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Since the temperature at fixed points (x p , y p ) was measured experimentally along the height of the 
enclosure (Ref. 1), Figure 12 provides the time history of the temperature field for various points along 
the enclosure height, contrasting low (Q = 0.01 W) and moderate (Q =1.0 W) heat absoiption. The 
experimental temperature measurements (Ref. 1 ) reported were averaged values over a plane consisting of 
five points at each y* location — near 0.1, 0.35, 0.65, and 0.90 — over a time span of t = 600 s. Figure 12 

shows boundary layer locations at y* p = 0.001 and 0.999 near the bottom and top wall boundaries to 
illustrate the instantaneous response of the fluid so that it can be compared with experimental trends. For 
Ra = 10 10 , in response to heat absoiption at the bottom boundary ( y* = 0.001), there is a shaip rise in 

temperature over a time scale of t = 100 s, and the onset of oscillations follows for t > 1 00 s. As the 
convective instability propagates toward the top of the enclosure, there is a delay time x D for the onset of 
convective instability: for y* p = 0.1, x z > = 150 s; whereas for y* p = 0.9, x D = 750 s. 

Similar to results for the bottom boundary, heat absoiption at the top boundary increases the local 
temperature; however, owing to the stably stratified layer on top (y* = 0.999), the temperature continues 

to rise steadily past t = 100 s up to t = 800 s, beyond which the onset of convective instability begins. 
Based on the delay time of x D = 800 s for convective instability to propagate the entire length of the 
enclosure, a heat load of Q = 0.01 W is not sufficient. For Ra = 10 10 (Q = 0.01 W), the delay time of 
Xd = 750 s at y * = 0.9 is too long for the propagation of convective instability in comparison to the 

experimental (Ref. 1) measurement time band of t = 600 s. This implies that the instability would not 
have reached the top of the enclosure. Thus, a background heat load of Q = 0.01 W is too small. Flowever, 
a decade increase to Ra = 10 11 (Q = 0.1 W) indicates a delay time of x D = 400 s at y * = 0.9, which falls 

within the experimentally measured time band of t = 600 s. Thus, Ra = 10 11 (Q = 0. 1 W) can potentially 
serve as a lower bound for heat absorption from the ambient laboratory environment. Flowever, when 
temperature magnitudes are compared with experimental values for Ra = 10 11 , they fall below the 
experimental measurements (Ref. 1). Another decade increase in heat load to Q = 1.0 W, Ra = 10 12 , 
shows a gradual approach to reported experimental values of average temperature and serves as a lower 
bound for heat absoiption from the ambient environment. 

The dynamics of the temperature field at fixed points (x p , y p ) spanning the height of the enclosure in 
Figure 12, shows a decrease in the penetrative convection response time as the Rayleigh number 
increases. In comparison to Ra = 1 0 10 , the delay response time of penetrative convective instability is 
reduced for Ra = 10 12 as the height increases. The response time ranges from x D = 25 to 155 s, corre- 
sponding to y* p =0.1 and 0.9, respectively. As the probing distance (y * ) increases from the bottom, the 

signal varies from a stationary time series about a constant mean value (y* p = 0.001) to a nonstationary 
time series about a mean value that increases with time (y* p = 0. 1). This trend is maintained as the 
enclosure height is traversed to the top boundary (y* p = 0.9), which indicates intermittency. In addition, 
the convective instability penetrates the top boundary layer at y* = 0.999, showing intermittency through 

periodic bursts of intensive convective motion. Owing to the lag in response time x D associated with 
locations from y* p =0.1 to 0.9 {x D = 25 to 155 s), a comparison with experiments over the short time 

scale of t = 60 s is inadequate for y* p > 0.3 since the experimental measurements indicate nearly instan- 
taneous rise in temperature as shown at y* p = 0.001 and 0.999. Flowever, for a time band of t > 200 s, 

adequate comparisons can be made with experimental data. 

The computational model accounts for the response of the temperature field to a step input of heat 
absoiption from the boundaries and represents a clearly posed initial value problem. Flowever, there exists 
uncertainty associated with the initial condition of an experimental system. Since heat is absorbed 
instantaneously, the origin of time ( t = 0) may be uncertain, considering the complexity of temperature 
measurements and the establishment of the initial saturation condition for liquid nitrogen. Considering the 
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Figure 12. — Comparison of low (heat load, Q = 0.01 W, and heat flux ratio at the sidewall, f s = 2.87x10~ 6 ) and 
moderate (Q = 1.0 W and f s = 2.87x1 0 -8 ) heat absorption from the ambient environment on penetrative 
convective instability along the height of the enclosure, as well as the typical dynamical state for Rayleigh 
number, Ra = 10 12 ; aspect ratio, Ar = 3.24 (height and length of enclosure, H = 69.1 cm and L =19.1 cm); 
Prandtl number, Pr = 2.21\ and heat flux ratios at bottom and top of enclosure, fb = 1 and f t = 1 , 
respectively. T, temperature; f, time; g", heat flux; k, constant. Results are shown at different points inside 
the enclosure (x*,y*). 
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complexity of an experiment, we conjecture that there could be uncertainty associated with the delay time 
x D for penetrative convection that cannot be captured in an experimental system. Nevertheless, as is 
shown in Section 4.6, the temperature magnitude over the time scale of the experiment, t = 600 s, 
indicates that the absoiption of a heat load of Q = 1.0 W from the background environment can serve as a 
lower bound for a potential scenario of the actual experiment. These results imply that the temperature 
signals from an experimental system, prior to averaging, are likely to be nonstationary, with a mean value 
that varies with time; for the experimental system (Refs. 1 and 2) considered, only averaged values of 
temperature were reported. A nonstationary time series requires proper treatment of the mean value 
(Ref. 12). Without proper averaging of a nonstationary time series, uncertainty can be introduced in the 
average value of temperatures. 

The time series shown for Ra = 10 12 at y* p = 0.9 is quantified using the amplitude spectrum of the 
u component of the velocity field, defined as 


U(J)~ 

1 R 


Tr 

J u{t)e~ 2m f l 

o 


(24) 


In relation to Figure 12 showing the dynamical state for Ra = 10 12 , Equation (24) shows the square root of 
the local kinetic energy per unit time, or the square root of power ( PJJ )) along the ordinate axis; the 
abscissa represents the spectrum of the frequency of the time signal; and T R denotes the interval of the 
time series. A convolution relationship, U(f) = W/,(f) * U if), using a Hanning spectral window W n (f) is 
used to obtain a smooth estimate U(f). The frequency decomposition of U(f) in Figure 12 indicates a 
long-period event that has a low-frequency peak of 0.005 Hz and spans a broadband of 0.5 Hz. The 
broadband response falls below the critical Nyquist frequency of f c = 3.5 Hz, indicating that the effective 
bandwidth of the signal is free from aliasing. Note that the power spectrum, P u (f) = U(f ) x U*(f)/2n, is 
the product of U(f) with its conjugate U*(f), which in terms of Equation (24) is the square of the 
magnitude of the Fourier-transformed velocity component u(t). The power spectrum gives results that are 
similar to those for the amplitude spectrum. However, the amplitude spectrum was found to have the 
advantage of showing the higher harmonics of the fundamental frequency on a linear scale. A 
pseudophase space trajectory, obtained by lagging the time signal of u(t + k) by a constant value k, shows 
a dense trajectory. The combination of a broadband power spectrum and a dense pseudophase space 
trajectory indicates a chaotic state. We conjecture that this chaotic dynamical state exists for the entire 
range of Rayleigh numbers considered — even for the lowest heat absorption (Q = 0.001 W) at Ra = 10 9 — 
owing to the large size of the enclosure. 

Related works on liquid helium by Behringer (Ref. 13) aimed at characterizing the onset of low- 
dimensional turbulence or chaos show that the trend of the critical Rayleigh number as a function of 
aspect ratio merges with the prediction of Charlson and Sani (Refs. 14 and 15) for long vertical walls, 
which shows an increase of the critical Rayleigh number for the onset of instability. In contrast, 
enclosures with large horizontal extent show that turbulence occurs near the critical Rayleigh number. 
These observations parallel our conjecture of chaos at the lowest Rayleigh number since our enclosure 
has relatively long vertical and horizontal length scales. 


4.6 Comparison of Model With Experiment 

Figure 1 3 compares the mean fluctuating temperature T at fixed points (x p , y p ) along the height of the 
enclosure with experimental data (Experiments (Exp.) — no screen, no heater) after Bolshinskiy et al. 

(Ref. 1). The model uses the boundary condition of heat load from Q= 1.0 to 10 W (Ra = 10 12 to 10 13 ) 
due to heat absoiption from the ambient laboratory environment. The experimental data (o) correspond to 
average values of temperatures <T> for five points over a plane at each height (y) location. We compare 
the experimental data with the mean temperature T at a fixed location in order to show unfiltered 
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predicted temperatures without plane averaging; note the unfiltered temperatures in Figures 1 1 and 12 
reveal inherent fluctuations that may not be captured by averaging. The experimental data points are in 
close proximity to the computational model, and for comparison, the figures show data for a long time, 
t = 600 s at Ra = 1 0 12 , and a short time, t = 50 s for Ra = 10 13 . Note that the delay time x D = 15 s for 
Ra = 10 13 is very short for high heat flux absoiption and falls within the response time of the experimental 
data. For Ra = 10 12 , a heat load of ()= I W yields temperature magnitudes that are lower than those 
measured in the bulk liquid; however, in the boundary layer near the top wall, the mean temperature, 

T = 19.5K, approaches the bulk experimental temperature, <T> = 80.8 K, of the liquid in the middle 
region. 

The model predicts the trend shown by the data: that is, a nearly uniform temperature in the central 
region of the enclosure, a nearly adverse temperature gradient near the bottom of the enclosure due to 
local convective instability, and a stable temperature gradient near the top region due to stable stratifi- 
cation. For an ambient heat load of Q = 10 W (Ra = 10 13 ) in the neighborhood of the maximum value 
estimated from thermal radiation due to the background laboratory environment (T x = 293 K), the 
saturation temperature of liquid nitrogen, T A = 77 K, nearly matches the experimental data for short times. 
In this case, the temperature in the boundary layer near the bottom wall, T = 77.5 K, nearly matches the 
highest temperature measured experimentally near the top region, <T> = 77.6 K. In contrast, the tempera- 
ture in the boundary layer near the top wall, T = 83.5 K , exceeds the highest measured temperature in the 
top region, <T> = 82 K, for the long duration time of t = 600 s. Given the potential uncertainty associated 
with establishing the proper initial condition that can occur for experimental systems, a minimum-to- 
maximum heat load range due to a background environment of Q = 1 to 10 W (Ra = 10 12 to 10 13 ) can be 
representative of actual experimental conditions. 



Figure 13. — Comparison of mean fluctuating temperature T at fixed points (x p , y p ) in the flow field with experimental 
data (Exp. — no screen and no heater) (Ref. 1 ) for increasing Rayleigh numbers, Ra. (a) Time, t = 600 s. (b) Time, 
t = 50 s. 
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5.0 Summary and Conclusion 

A model problem is addressed on the effect of heat absorption from the background laboratory 
environment on the temperature rise of liquid nitrogen at saturation conditions (p^, T ,) inside a large 
enclosure (H = 61.9 cm and L= 19.1 cm). The model contributes toward an understanding of the flow 
field dynamics for an experimental system aimed at investigating heat entrapment effects in liquid 
acquisition devices. For estimated background heat load conditions, we consider a range from Q = 0.001 
to 10W, corresponding to Rayleigh numbers fromRn = 10 9 to 10 13 . The estimated heat flux of the top 
and bottom boundaries ( q" = ql ) ranges from 3.49x10' to 3.49xl0 5 erg/cm 2 -s, and a conservative 
estimate of q" = 0.001 erg/cm 2 -s is used for the sidewalls owing to the use of insulation at the side 
boundary in the experiments. In response to heat load, the flow field dynamics for Ra = 10 9 to 10 13 show 
self-similarity. The flow field is driven by local convective instability, against a free convection base 
flow, generated at the bottom of the enclosure. The flow field shows small-scale thermals for the short 
time scale and a traveling-wave -type motion of convective modes for the long time scale. 

Finite heat leakage at the sidewalls generates the initial free-convection base flow, which consists of 
two counterrotating cells that become unstable against convective instability caused by an adverse 
temperature gradient at the bottom boundary, while the top boundary forms a stably stratified configura- 
tion. A number of modes arise, depending on the Rayleigh number, which self-organize into two modes 
owing to amalgamation prior to the propagation of the convective instability along the height of the 
enclosure. A number of symmetric vertical modes arise; these modes become asymmetric prior to 
reaching the top of the enclosure depending on the Rayleigh number. A traveling-wave-type motion of 
convective modes follows that generates aperiodic oscillations of the flow and temperature fields and is 
representative of the asymptotic dynamics. 

The time history of the temperature and velocity fields shows a characteristic delay time x D for the 
penetration of the convective instability along the height of the enclosure due to heat load absorption. The 
estimated heat loads, which serve as lower and upper bounds, owing to the background environment 
range from Q = 1 to 10 W, and show agreement with experimental trends (Ra = 10 12 to 10 13 ). For heat 
loads from Q = 1 to 1 0 W, the time signal is nonstationary over a time band and spatial location range, 
with a mean that increases with time. Thus, cautionary measures should be taken for averaging of 
nonstationary temperature measurements. A chaotic dynamical state (stochastic time signal) seems to be 
the norm for the range of Rayleigh numbers considered, Ra = 10 9 to 10 13 , owing to the large size of the 
enclosure. 

The computational model sheds light on temperature measurements of saturated liquid nitrogen inside 
an enclosure subjected to an unknown heat load from a background laboratory environment for applica- 
tion to liquid acquisition devices. To understand the effect of the applied heat flux at the bottom of the 
enclosure on the temperature rise of liquid nitrogen, it is necessary to establish the base condition of the 
experiment. This base condition, which was characterized experimentally by measurements of the 
temperature rise due to the background environment (Exp. — no screen, no heater) can be used to establish 
the boundary conditions. The computational model shows that the boundary conditions for the experiment 
lie between Q = 1 and 10 W (Ra = 10 12 and 10 13 ). These established boundary conditions can be used as 
input for the top heat flux q" along with the estimated side heat flux q" in order to analyze the next phase 
of the experiment, that is, the effect of independent variation of the bottom heat flux q” h . This implies 
that a bottom heat-load variation from a minimum of Q = 1 0 W to the maximum experimental input of 
Q = 1 04 W can be investigated. 
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aspect ratio 
heat flux ratio 
critical Nyquist frequency 
gravitational acceleration 
body force in y direction 
height of enclosure 
time lag constant 
horizontal length of enclosure 

coefficient representing gravitational level, where n = 1 represents the gravity on Earth 
power spectrum of u velocity component 
Prandtl number 
pressure 

heat load from ambient 
imposed heat flux 
Rayleigh number 
temperature 

mean fluctuating temperature 
temperature of saturated liquid nitrogen 
interval of the time series 
ambient temperature of laboratory 
time 

characteristic time for formation of thermals 

smooth estimate of the amplitude spectrum of the u component of the velocity field (after 
Eq. (24)) 

amplitude spectrum of the u component of the velocity field (Eq. (24)) 

horizontal component of velocity 
velocity magnitude 

velocity field 

vertical component of velocity 
Hanning spectral window 
fixed horizontal point 
fixed vertical point 
thermal diffusivity 
thermal expansion coefficient 
boundaiy of enclosure 
thermal boundaiy layer thickness 
thermal conductivity 
kinematic viscosity 

single vorticity component in the z direction (Eq. (11)) 


N ASA/TM— 20 14-21 7442 


23 



p density 

p A density of the saturated liquid nitrogen 

x d delay time 

\| / stream function 

V for all 

V gradient operator 

A parametric set (Eq. (22)) 

A' reduced parametric set (Eq. (23)) 

i magnitude 

Subscripts: 

b bottom 

p pressure (after Eq. (3)) 

s sidewall 

t top 

Superscripts: 

* dimensionless quantity 
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